coeftest.cluster <- function(data, fm, cluster1=NULL, cluster2=NULL, ret="test") {
  
  library(sandwich)
  library(lmtest)
  
  data <- as.data.frame(data)
  
  # Return White (1980) standard errors if no cluster
  # variable is provided
  if (is.null(cluster1)) {
    if (ret=="cov") {
      return(vcovHC(fm, type = "HC0"))
    } else {
      return(coeftest(fm, vcov = vcovHC(fm, type = "HC0")))
    }
  }
  
  # Calculation shared by covariance estimates
  est.fun <- estfun(fm)
  
  # Need to identify observations used in the regression (i.e.,
  # non-missing) values, as the cluster vectors come from the full
  # data set and may not be in the regression model.
  inc.obs <- !is.na(est.fun[,1])
  est.fun <- est.fun[inc.obs,]
  
  # Shared data for degrees-of-freedom corrections
  N  <- dim(fm$model)[1]
  NROW <- NROW(est.fun)
  K  <- fm$rank
  
  # Calculate the sandwich covariance estimate
  cov <- function(cluster) {
    cluster <- factor(cluster, exclude=NULL)
    
    # Calculate the "meat" of the sandwich estimators
    u <- apply(est.fun, 2, function(x) tapply(x, cluster, sum))
    meat <- crossprod(u)/N
    
    # Calculations for degrees-of-freedom corrections, followed
    # by calculation of the variance-covariance estimate.
    # NOTE: NROW/N is a kluge to address the fact that sandwich
    # uses the wrong number of rows (includes rows omitted from
    # the regression).
    M <- length(levels(cluster))
    dfc <- M/(M-1) * (N-1)/(N-K)
    
    #print (sandwich(fm, meat=meat))
    return(dfc * NROW/N * sandwich(fm, meat=meat))
  }
  
  # Calculate the covariance matrix estimate for the first cluster.
  cluster1 <- data[inc.obs, cluster1]
  cov1  <- cov(cluster1)
  # print(cov1)
  
  if (is.null(cluster2)) {
    # If only one cluster supplied, return single cluster
    # results
    if (ret=="cov") {
      return(cov1)
    } else {
      return(coeftest(fm, cov1))
    }
  } else {
    # Otherwise do the calculations for the second cluster
    # and the "intersection" cluster.
    cluster2 <- data[inc.obs, cluster2]
    cluster12 <- paste(cluster1, cluster2, sep="")
    
    # Calculate the covariance matrices for cluster2, the "intersection"
    # cluster, then then put all the pieces together.
    cov2   <- cov(cluster2)
    cov12  <- cov(cluster12)
    covMCL <- (cov1 + cov2 - cov12)
    
    # Return the output of coeftest using two-way cluster-robust
    # standard errors.
    # print(ret)
    if (ret=="cov") {
      return(covMCL)
    } else {
      return(coeftest(fm, covMCL))
    }
  }
}

# Following based on suggestion from
# https://stat.ethz.ch/pipermail/r-help/2011-January/264777.html
# provided by Achim Zeileis.
summary.cluster <- function(obj, data, cluster1, cluster2=NULL, alpha=0.05) {
  
  require(memisc)
  
  # Get original summary
  s <- getSummary(obj, alpha=alpha)
  
  ## replace Wald tests of coefficients
  s$coef[,1:4] <- coeftest.cluster(data, obj, cluster1, cluster2)
  
  ## replace confidence intervals
  crit <- qt(alpha/2, obj$df.residual)
  s$coef[,5] <- s$coef[,1] + crit * s$coef[,2]
  s$coef[,6] <- s$coef[,1] - crit * s$coef[,2]
  
  # Note that some components of s$sumsstat will be inconsistent with
  # the clustered calculations
  
  return(s)
}